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We  have  obtained  realistic  fluid  detonation-wave  profiles.  1 

This  work  is  summarized  in  Pages  2-18. 

I 

We  have  studied  the  structure  of  uni  axially  and  hydro¬ 
statically-compressed  solids  and  the  transfer  of  energy  among 

> 

translational  and  internal  molecular  modes.  This  work  is  described 
in  Pages  19-26. 

We  have  developed  novel  computational  methods  for  simulating 
nonequilibrium  processes  using  Gauss'  Principle  of  Least  Constraint. 

This  work  is  described  in  Pages  27-38. 
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SUMMARY 


(1)  TECHNICAL  PROBLEM 

We  are  investigating  mathematical  models  of  detonation  in  dense  fluids  and 
in  solids.  We  wish  to  obtain  a  fundamental  understanding  of  the  interaction  of 
pressure  waves,  or  shockwaves,  with  high  explosives. 

(2)  GENERAL  METHODOLOGY 

We  are  solving  the  steady-state  dynamical  equations  of  continuun  mechan¬ 
ics,  including  viscosity,  conductivity,  and  chemical  reactions,  to  simulate 
detonation  in  fluids  capable  of  undergoing  fast  exothermic  reactions.  We  are 
solving  the  coupled  ordinary  differential  equations  of  motion  of  molecular 
dynamics,  for  solids,  to  simulate  the  initiation  process  preceding  detonation. 
We  are  examining  simple  models  designed  to  describe  the  results  of  these  macro¬ 
scopic  and  microscopic  approaches. 

(3)  TECHNICAL  RESULTS 

We  have  obtained  detonation  wave  profiles,  using  a  realistic  dense-fluid 
equation  of  state.  We  have  compared  these  profiles  to  the  predictions  of  the 
simplified  Zeldovich-Von  Neumann-Doering  model. 

We  have  developed  methods  for  simulating  the  rapid  uniaxial  compression  of 
solids  suited  to  molecular  dynamics  simulation. 

We  have  formulated  the  interatomic  distribution  function  in  solids,  for 
comparison  with  atomistic  simulations  and  use  in  kinetic  reaction  models. 

(4)  IMPLICATION  FOR  FURTHER  RESEARCH 

The  relatively  long  vibrational  relaxation  time^  occurring  in  shocked 
diatomic  molecules  with  realistic  values  of  viscosity  and  thermal  conductivity 
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suggest  that  the  nonequil Ibrlum  distribution  of  energy  is  more  important  than 
the  interaction  of  ordinary  transport  properties  with  chemical  reactions.  We 
therefore  intend  to  concentrate,  in  the  next  two  years,  on  reaction  initiation 
and  intramolecular  energy  transfers  in  reacting  molecules  in  the  solid  phase. 


I.  INTRODUCTION 


Shock  and  detonation  waves  involve  atomic-scale  processes  which  are 
Imperfectly  understood  and  hard  to  measure.  For  this  reason  we  have  embarked 
on  a  program,  using  modern  computational  methods,  of  studying  the  processes 
through  which  mechanical  shockwave  energy  stimulates  the  release  of  chemical 
energy  within  a  detonation  wave. 

Our  work  entails  the  correlation  of  two  points  of  view:  continuum 
mechanics  and  atomic  mechanics.  We  seek  to  understand  the  limitations  of  the 
continuun  approach  in  dealing  with  atomic-scale  problems  and  to  develop  methods 
for  describing  these  small-scale  processes  in  a  realistic  way. 

Considerable  information  on  the  equilibrium  and  transport  properties  of 
simple  systems  is  now  avail ab 1 e ( 1 ,2) .  Computer  simulations  of  dense  fluids, 
coupled  with  thermodynamic  perturbation  theories(3)  based  on  simple  hard  and 
soft  sphere  models,  make  it  possible  to  describe  the  pressure-vol ume- 
temperature  equation  of  state  accurately,  for  sufficiently  simple  materials. 

For  some  "realistic"  potentials,  such  as  the  Lennard-Jones  6-12  interatomic 
potential,  comprehensive  computer  simulation  studies  have  established  the 
density  and  temperature  dependence  of  the  viscosities  and  the  thermal  conducti¬ 
vity  as  well(2).  Much  less  is  known  about  the  contribution  of  intramolecular 
degrees  of  freedom  to  the  transport  properties  and  equation  of  state.  Studies 
on  molecules  like  methane  and  benzene  indicate  that  the  computer  simulation 
approach  is  now  feasible,  even  for  molecules  containing  several  atoms(4). 

A  numerical  study,  comparing  the  continuum  and  atomistic  versions  of  a 
very  strong  dense-fluid  shockwave  made  use  of  both  the  equilibrium  and  non- 
equilibriun  constitutive  properties  of  the  6-12  potential,  and  indicated  that 
useful  conclusions  could  be  drawn  by  comparing  the  two  approaches(5) .  A 
natural  extension  of  that  work  is  to  consider  the  possibility  of  coupling 


chemical  reactions  with  the  dissipative  processes  of  viscosity  and  thermal 
conductivity  which  appear  in  shockwaves.  That  coupled  combination  is  a  detona¬ 
tion  wave. 

The  simplest  prototype  condensed-phase  high  explosive  is  probably  nitric 
oxide,  NO.  This  molecule  forms  molecular  nitrogen  and  oxygen  in  an  exothermic 
reaction  generating  approximately  100  kilobars  pressure(6).  We  set  out  with 
the  goal  of  simulating  the  detonation  of  liquid  NO  using  both  the  continuum  and 
atomistic  approaches.  The  continuum  calculations  are  complete.  Using  an 
equation  of  state  which  describes  the  Lennard-Jones  potential  throughout  the 
fluid  regions  of  the  phase  diagram,  and  transport  coefficients  from  the  Enskog 
theory,  we  have  generated  detonation  profiles  as  functions  of  the  energy 
release  and  activation  energy  of  the  dense-fluid  detonation.  These  calcula¬ 
tions,  including  viscosities  and  thermal  conductivity,  are  compared  to  the 
Zeldovich-Von  Neumann-Doering  model  in  Section  II. 

Shortly  after  this  work  was  begun,  a  research  program  was  initiated  at  Los 
Alamos,  designed  to  measure  and  to  model  the  detonation  of  nitric  oxide.  This 
program  includes  experimental,  quantum  mechanical,  and  dynamical  groups  working 
toward  a  common  goal.  Related  Los  Alamos  research,  carried  out  by  Brad  Hoi i an, 
Galen  Straub,  and  their  co-workers,  has  established  that  the  equilibration  time 
between  vibrational  and  translational  temperatures  in  dense  diatomic  fluids 
(nitrogen)  is  relatively  long(7),  just  as  in  the  gas  phase.  These  long  equili¬ 
bration  times,  of  order  100  vibration  times,  suggest  that  the  complete  combus¬ 
tion  process  is  impractical ly  long  for  an  atomistic  simulation.  For  this 
reason  we  have  concentrated  our  more  recent  work  on  solid-phase  initiation  and 
mode  equilibrium.  We  are  developing  solid-phase  calculations  in  which 
chemically  active  molecules  can  be  studied  under  uniaxial  and  homogeneous 


adiabatic  compression.  That  work  is  described  in  Sections  III  and  IV.  We 
expect  to  devote  the  next  two  years'  effort  to  enhancing  our  understanding  of 
solid-phase  detonation  initiation  from  the  atomistic  mechanical  point  of  view 


II.  DENSE-FLUID  DETONATION  WAVE  STRUCTURE 

The  classical  steady  one-dimensional  ZND  (for  Zeldovich,  Von  Neumann,  and 
Doering)  model  of  detonation-wave  structure  can  be  summarized  with  the  help  of 
Figure  1.  It  is  assumed  that  the  reactants  can  assume  the  equilibrium  unreact¬ 
ed  thermodynamic  states  shown  on  the  lower  Hugoniot  curve  as  the  result  of 
shock  compression.  Each  of  these  states  satisfies  the  condition  that  its  mass, 
momentum,  and  energy  fluxes  agree  with  those  characterizing  the  initial  state 
(state  0).  This  unreacted  Hugoniot  is  hypothetical  because  chemical  reactions 
will  gradually  transform  the  reactants,  at  least  for  strong  enough  shockwaves, 
to  hot  products. 

The  fully  reacted  products  can  also  be  described  by  a  Hugoniot  curve 
satisfying  the  conservation  laws.  This  is  the  upper  Hugoniot.  Its  position, 
relative  to  the  lower  Hugoniot,  depends  upon  the  heat  of  reaction  Q.  The 
Zeldovich-Von  Neumann-Doering  detonation  wave  model  pictures  the  detonation 
process  as  two  consecutive  steps:  first,  the  reactants  are  shock-compressed  to 
the  (Von  Neumann  spike)  state  S.  Next,  these  hot  reactants  slowly  undergo 
chemical  reactions,  proceeding  to  the  final  (Chapman-Jouguet)  state  CJ.  During 
the  reaction  phase  the  state  point  follows  the  Rayleigh  line,  satisfying  the 
constraints  of  constant  mass  and  momentum  flux  necessary  in  a  steady  wave. 
Neither  the  transport  coefficients  nor  the  reaction  rate  enter  into  Figure  1. 
These  state  functions  determine  the  space  and  time  scales  associated  with  the 
compression  and  reaction  phases  of  the  detonation  wave. 

The  ZND  model  is  an  oversimplification.  With  the  theory  of  liquids  suf¬ 
ficiently  advanced  that  quantitative  calculations,  including  nonideal  effects, 
can  be  made,  there  is  no  real  need  for  an  inaccurate  treatment. 

Molecular  dynamics  simulations  have  shown  that  strong  dense-fluid  shock- 
waves  can  be  described  fairly  well  by  the  Navi er-Stokes  equation(5).  Figure  2 
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Indicates  the  relatively  good  agreement  between  the  Navi er-Stokes  continuum 
calculations  and  the  atomistic  simulations.  The  principle  disagreements  are 
these:  The  wave  is  about  30%  broader  than  the  continuum  predictions  and  the 
velocity  distribution  is  considerably  more  complex  than  is  the  equilibrium 
Maxwel 1-Boltzmann  distribution.  The  alternative  Mott- Smith  approach(8),  in 
which  the  velocity  distribution  is  a  weighted  average  of  the  reactant  and 
product  Maxwellian  distributions,  grossly  overestimates  the  nonequilibrium 
character  of  the  wave. 

We  initially  planned  to  carry  out  a  similar  comparison  for  nitric  oxide,  a 
simple  liquid  which  undergoes  the  reaction 

2N0  =  N2  +  02.  (1) 

Exploratory  calculations  fitted  to  the  kinetics  of  the  reaction  showed  that  the 
reaction  zone  extends  for  thousands  of  angstroms,  far  too  long  for  molecular 
dynamics  simulation.  The  liquid  detonation  wave  profile  that  we  display  in 
Figure  3,  together  with  the  ZNO  approximation  to  that  profile,  was  calculat¬ 
ed^)  with  an  activation  energy  approximately  half  that  of  nitric  oxide.  This 
substantially  reduced  the  spatial  extent  of  the  detonation  wave.  Under  this 
artificial  assumption  there  is  a  noticeable  quantitative  discrepancy  between 
the  simple  two-step  model  and  the  solution  of  the  coupled  equations.  There  is 
no  such  dramatic  difference  for  nitric  oxide,  at  least  for  a  one-dimensional 
wave. 


The  instability  of  one-dimensional  low-density  detonation  waves  is  well- 
established  experimental ly( 10) .  Because  this  effect  is  at  least  as  important 
as  that  of  the  hydrodynamic  transport  coefficients  we  see  that  further  study  of 
the  one-dimensional  detonation  wave  for  nitric  oxide,  is  not  feasible,  using 
molecular  dynamics. 


Vibrational  energy  transfer  Is  very  slow  In  diatomic  molecules  because  the 
single  vibrational  mode  is  well-separated  from  the  "lattice  modes."  The  situa¬ 
tion  is  less  clear  cut  for  polyatomic  molecules  with  many  degrees  of  freedom. 
Can  the  solid-phase  initiation  problem  be  understood  using  equilibrium  chemical 
kinetics  and  equilibrium  energy  surfaces?  Are  nonequilibrium  velocity  distri¬ 
butions  and  intramolecular  energy  transfer  significant?  The  answers  to  these 
questions  are  unknown  despite  considerable  interest  and  debate(ll).  We  are 
concentrating  our  efforts  on  enhanced  understanding  of  solid-phase  Initiation, 
using  the  methods  outlined  in  the  following  sections. 


■■  p  m  ■  y  >  j*  '."B'j  ■  »'i 


T»1 


III.  COMPUTER  SIMULATIONS  AT  HIGH  RATES  OF  STRAIN 

To  explore  the  rapid  deformation  associated  with  shockwaves  in  condensed 
fluids  and  solids  we  have  developed  special  methods  for  solving  the  equations  of 
motion  of  atoms  in  a  rapidly-deforming  fluid  or  sol id(12) .  These  methods  were 
originally  developed  and  applied  to  the  simulation  of  viscous  flows(13)  and  were 
later  applied  to  solid-phase  pi asticity( 14) .  In  either  case  the  macroscopic 
flow  is  described  by  a  given  strain-rate  tensor  Vu,  where  u  is  the  hydrodynamic 
flow  velocity.  To  incorporate  the  flow  velocity  into  the  equations  of  motion,  a 
perturbation,  the  tensor  product  of  Vu  and  Doll's  tensor  (the  dyadic  composed 
the  particles'  summed  coordinate-momentum  product  qp)  is  added  to  the 
Hamiltonian  function.  Equations  of  motion  result  which  include  contributions 
proportional  to  the  strain  rate: 
q  =  (p/m)  +  q-Vu  ; 

p  =  F-  Vu-p;  (2) 

E  =  -  VP:VU  . 

Equations  (2)  describe  adiabatic  homogeneous  deformation.  This  scheme  has 
been  used  to  make  quantitative  estimates  of  the  shear  and  bulk  viscosities  for 
dense  fluids. 

Inhomogeneous  compressions  must  be  treated  differently.  The  system  is 
sheared  or  compressed  by  moving  periodic  images(5).  Either  method  could  be 
applied  to  chemically-reacting  solids.  We  expect  to  explore  the  homogeneous 
uniaxial  compression  of  two  systems:  The  simple  reaction  model  described  in  the 
next  section  and  a  "realistic"  three-dimensional  polyatomic  model  of  a  detonat¬ 
ing  solid.  This  latter  application  should  allow  us  to  map  out  the  path  followed 
by  energy  during  the  equi 1 ibration  process. 
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IV.  ATOMIC  DISTRIBUTION  FUNCTIONS  IN  SOLIDS 


We  are  currently  investigating  a  simple  model  for  reacting  atoms:  Two 
atoms  sufficiently  close  together  react,  releasing  an  energy  Q.  Equilibrium 
simulations  (with  Q=0)  have  been  carried  out  to  determine  the  density,  strain, 
and  temperature  dependence  of  the  "reaction  rate."  Figure  4  indicates  the 
substantial  dependence  on  strain  biaxiality.  Pastine's  model(15),  which 
contains  an  adjustable  collision  frequency  w,  can  be  used  to  fit  the  uniaxial 
compression  results. 

A  more  fundamental  approach  incorporates  the  distribution  of  neighboring 
atoms.  In  the  quasiharmonic  case  this  can  be  done  exactly.  We  have  developed 
a  numerical  method  for  the  partial  integration  of  the  Boltzmann  factor, 
exp(-4|/kT) ,  over  all  but  two  particle  coordinates,  leading  to  an  atom-atom 
distribution  function.  This  function  is  a  product  of  Gaussians. 

Because  t  calculation  is  somewhat  intricate  it  is  important  to  have  a 
check  available.  The  close-packed  lattice,  with  nearest-neighbor  interactions, 
is  ideal  for  this  purpose.  Because  the  potential  energy  of  a  nearest-neighbor 
pair  in  a  D-dimensional  close-packed  crystal  Is  kT/(D+l),  the  mean-squared 
deviation  in  the  nearest-neighbor  spacing  should  be  2kT/(D+l)<,  where  <  is  the 
nearest-neighbor  force  constant.  We  have  verified  that  this  result  is  repro¬ 
duced  exactly  by  our  numerical  approach,  and  find  that  the  convergence  of 
higher-order  and  transverse  correlations  is  relatively  rapid  and  smooth  as  the 
number  of  particles  in  the  periodic  cell  is  increased.  In  the  two-dimensional 
case,  for  example,  the  ratios  of  the  nearest-neighbor  mean-squared  longitudinal 
to  transverse  fluctuations  for  9,  16,  25,  and  36-atom  crystals  are  0.667, 

0.688,  0.695,  and  0.696(16). 

We  expect  to  apply  this  same  technique  to  the  atom-atom  distribution 
function  for  polyatomic  crystals  undergoing  uniaxial  compression.  This  analog 
of  Pastine's  model’ should  prove  useful  in  estimating  the  dependence  of  reaction 
rate  on  the  macroscopic  strain  rate.  _n_ 
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FIGURES: 

1.  Diatomic  reactant  and  product  Hugoniots  calculated  with  the  Lennard-Jones 
equation  of  state.  The  well  depth  and  collision  diameter  are  e  and  o, 
respectively.  The  initial  volume  Is  0.8442Na3.  The  initial  temperature 
is  0.722e/k. 

2.  Shockwave  profile  for  a  strong  fluid  shock  carrying  the  triple-point  argon 
to  a  final  temperature  of  about  one  electron  volt.  The  smooth  curves  are 
the  Navier-Stokes  solution.  The  points  correspond  to  data  from  nonequili¬ 
brium  molecular  dynamics. 

3.  Detonation  profile  calculated  using  (a)  the  hydrodynamic  equations,  includ¬ 
ing  viscosity  and  heat  conduction,  and  (b)  using  the  ZND  model.  The  activa- 
tation  energy  is  artifically  reduced  to  enhance  the  difference  between  the 
two  approaches.  The  heat  of  reaction  is  100c  and  the  activation  energy  is 
70e. 

4.  Reaction  rate  as  a  function  of  density  for  a  simple  binary-collision  model. 
Both  uniaxial  and  hydrostatic  compressions  are  shown  at  a  temperature  about 
half  the  melting  temperature.  Pastine's  model  for  the  rate  is  shown  as  the 
full  curve. 
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SUMMARY 


(1)  TECHNICAL  PROBLEM 

We  are  investigating  atomistic  mathematical  models  of  energy  transferred  in 
condensed  media.  We  wish  to  obtain  a  fundamental  understanding  of  the  inter¬ 
action  of  compressive  waves  with  intermolecul ar  and  intramolecular  energies, 
leading  to  the  initiation  of  detonation. 

(2)  GENERAL  METHODOLOGY 

We  analyze  quasistatic  compression  by  comparing  equilibrium  statistical- 
mechanical  calculations  with  the  results  of  equilibrium  Newtonian  molecular 
dynamics  simulations.  We  are  developing  and  solving  nonequilibrium,  equations  of 
motion  to  describe  the  rapid  compression  of  microscopic  polyatomic  systems.  We 
emphasize  the  flow  of  energy  among  the  molecular  degrees  of  freedom  excited  by 
the  compression  process. 

(3)  TECHNICAL  RESULTS 

We  have  determined  the  atomistic  pair  distribution  function  for  planar 
crystals  undergoing  both  uniaxial  and  hydrostatic  compression.  We  have  compared 
these  distributions  to  the  predictions  of  the  Pasti ne-Kamlet- Jacobs  model. 

We  have  simulated  bimolecular  collisions  of  tri atomic  and  hexatomic  planar 
molecules  (with  hexani trobenzene  in  mind).  We  analyze  the  normal-mode  vibrations 
before  and  after  collision  and  observe  the  transfer  of  energy  among  the  transla¬ 
tional,  rotational,  and  vibrational  modes. 

(4)  FURTHER  RESEARCH  IMPLICATIONS 


.  \  -V  -V.  .  .  -  ., 


I.  INTRODUCTION 

Shock  and  detonation  waves  involve  atomic-scale  processes  which  are  imper¬ 
fectly  understood  and  hard  to  measure.  The  systems  contain  too  many  atoms  for  an 
accurate  quantum-mechani cal  treatment.  Classical  calculations  can  deal  with 
systems  of  hundreds  or  thousands  of  interacting  atoms,  over  hundreds  or  thousands 
of  vibrational  periods. 

To  evaluate  the  utility  of  the  microscopic  approach  we  have  embarked  on  a 
program  to  develop  and  apply  modern  computational  methods  to  the  nonequilibrium 
dynamical  processes  responsible  for  explosive  initiation.  A  detailed  microscopic 
view  is  necessitated  by  the  violent  nonequilibrium  nature  of  shock  compression. 
Simulations  show  that  the  intrinsic  width  of  a  shock  front  is  only  a  few  atomic 
diameters. 1  Thus  the  usual  concept  of  equilibrium  temperature,  with  Arrhenius 
kinetics,  is  unlikely  to  provide  an  understanding  of  the  initiation  of  chemical 
reactions.  To  achieve  understanding  of  initiation  detailed  studies  of  rapid 
nonequilibrium  i ntermol ecul ar  and  intramolecular  energy  transfer  must  be  carried 
out  and  analyzed. 

Over  the  past  several  years  there  has  been  considerable  progress  in  treating 
the  deformation  of  simple  monatomic  fluids  and  solids.  The  techniques  developed 
have  been  tested  by  intercomparison  with  experimental  data  and  with  less- 
efficient  older  simulation  techniques."^  The  irreversible  heating  associated 
with  rapid  vi scons  and  plastic  deformation  has  been  avoided  by  developing  iso¬ 
thermal  equations  of  motion.^  These  equations  constrain  the  second  moment  of 
the  velocity  distribution,  <(v-<v>  >,  to  a  fixed  value.  Analogous  developments 
for  treating  diffusion  and  conduction  stem  from  this  same  foundation.^ 

These  new  methods  have  not  been  applied  to  flexible  polyatomic  molecules.  In 
order  to  carry  out  this  generalization  we  have  analyzed  both  quasistatic  and 
rapid  deformations  of  monatomic  crystals  (Section  II),  the  rotational  and  col  1 i - 
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II.  DEFORMATION  OF  MONATOMIC  SOLIDS 

In  a  series  of  studies  we  have  characterized  the  compression  of  simple 
fluids  and  solids.  We  recently  carried  out  a  detailed  study  of  equilibrium  pair 
distribution  functions  in  deformed  two-dimensional  crystals.®  Both  Hooke's-Law 
and  Lennard-Jones  forces  were  used.  The  pair  distributions,  obtained  by  integra¬ 
ting  the  quasi harmoni c  canonical  probability  density,  differ  from  approximate 
theoretical  treatments  in  two  ways.  First,  contributions  from  all  normal-mode 
vibrations  are  included.  Second,  reaction-rate  contributions  from  all  diatomic 
orientations  are  included.  For  these  reasons  the  reaction  rates  follow  neither 
the  simple  Gruneisen  density  dependence  nor  Arrhenius  kinetics.  Even  under 
quasistatic  conditions,  uniaxial  compression  generally  yields  significantly 
higher  collision  (or  reaction)  rates  than  does  the  corresponding  hydrostatic 
compression. 

Earlier  deformation  studies,  under  both  shockwave  and  homogeneous  condi¬ 
tions,  have  demonstrated  the  utility  of  incorporating  the  macroscopic  velocity 
gradientVu  directly  into  microscopic  equations  of  motions.  The  shockwave 
studies  emphasized  the  nonequi 1 i bri um  velocity  distribution  in  the  shock-front. 
Shockwave  structure,  bulk  and  shear  viscosities,  and  solid-phase  high-strai n-rate 
yield  stengths^  have  all  been  determined  in  this  way.  Adiabatic  deformation 
can  be  modeled  by  adding  an  extra  contribution  to  the  atomic  velocities, 
q  =  q-Vu,  and  a  correspondi ng  force  p  =  -7u» p.  The  two  modifications  lead  to 
the  thermodynamic  identity 

*E  =  -VP:7u, 

where  E  is  internal  energy  and  P  is  the  pressure  tensor.  This  adiabatic  deforma¬ 
tion  can  be  made  isothermal  instead  by  including  an  additional  collective  force 
-^p,  with'J  chosen  to  keep  <p^/m>  constant. 
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III.  ENERGY  TRANSFER  IN  POLYATOMIC  MOLECULES 


The  simplest  polyatomic  molecule  exhibiting  energy  transfer  is  a  planar 
triatomic  molecule,  with  three  vibrational  degrees  of  freedom.  In  this  case  the 
two  degenerate  vibrational  modes  are  coupled  by  Coriolis  forces.  Thus  rotation 
of  the  molecule,  even  very  slowly,  leads  to  a  geometric  coupling  of  two  vibra¬ 
tional  amplitudes. 

Regular  hexatomic  planar  molecules  have  nine  vibrational  degrees  of  free¬ 
dom,  with  four  pairs  of  degenerate  modes  in  addition  to  the  symmetric  breathing 
mode.  At  reasonable  rotational  velocities  the  coupling  among  the  four  pairs  is 
negligible.  Even  at  vanishing  velocities  Cordiolis  coupling  again  causes  the 
cycling  back  and  forth  between  degenerate  mode  amplitudes. 

We  have  simulated  bimolecular  collisions  of  both  triangular  and  hexagonal 
planar  molecules,  have  analyzed  their  normal-mode  vibrations  throughout  the 
collision  process,  and  have  observed  energy  transfer  among  all  the  translational, 
rotational,  and  vibrational  degrees  of  freedom.  Presently  we  are  characterizing 
the  energy  distribution  during  a  homogeneous  compression. 


IV.  THERMODYNAMICS  OF  RAPID  DEFORMATION 


The  virial  theorem  relates  the  pressure  tensor  to  interatomic  forces.  In 
the  familiar  monatomic  form®  each  Newton's  equation  of  motion  is  multiplied  by 
the  corresponding  particle  coordinate,  summed  and  averaged,  to  give 
Z^r-jF-jj  -  VP  =2'^mrr)i  -  (mrr)i. 
which  is  equivalent  to 

PV  =  NkTI  +ZirijFij 

where  N  is  the  number  of  atoms  in  the  volume  pairs.  For  polyatomic  molecules  it 
is  convenient  to  multiply  each  atom's  equation  of  motion  by  the  center  of  mass 
coordinate  for  its  molecule.  The  force  terms  in  each  molecule  cancel,  with  the 
result: 

PV  =  NkTI  +2>ijFij 

where  N  is  now  the  number  of  molecules  and  F. j  is  the  vector  sum  of  forces  of 
all  atoms  in  molecule  I  due  to  atoms  in  molecule  J. 

Now  consider  two  kinds  of  adiabatic  deformation.  In  "atomic"  deformation 
each  atom  undergoes  a  displacement  proportional  to  its  location 

r  =  r-Vu 


and  a  correspond! ng  force 

*p  =  -P'V  u. 

r 

The  adiabatic  first-law  identity,  E  =  -VP-Vu,  follows.  In  "molecular"  deforma¬ 
tion  each  atom  undergoes  a  di spl acment  proportional  to  the  center-of -mass 
coordinate  R: 

r  =  R  '  Vu 


as  well  as  a  corresponding  force 

P  =  -PC*)*  Vu. 

The  identify  E  r  -PV:Vu  again  follows,  provided  that  P  is  the  "molecular" 
pressure. 
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SUMMARY 


(1)  TECHNICAL  PROBLEM 

We  consider  atomistic  descriptions  of  matter,  far  from  equilibrium,  in  order 
to  understand  the  mechanisms  underlying  detonation  in  high  explosives.  We  focus 
on  solid  hexani trobenzene. 

(2)  GENERAL  METHODOLOGY 

We  develop  equations  of  motion  incorporating  both  thermal  and  mechanical  con¬ 
straints.  Adiabatic,  isothermal,  and  isoenergetic  processes  can  be  simulated  by 
solving  these  equations  numerically.  We  study  many  special  cases  involving  gases, 
liquids,  and  solids. 

(3)  TECHNICAL  RESULTS 

We  have  established  the  usefulness  and  the  limits  of  "Gauss'  Principle  of 
Least  Constraint"  in  classical  many-body  simulations.  We  have  applied  three  dif¬ 
ferent  types  of  feedback  control  to  atomistic  systems  in  studying  fluid  and  solid 
heat  flow.  We  have  considered  the  size-dependence  of  constitutive  properties.  We 
have  developed  a  simple  stress-free  model  for  crystalline  hexani trobenzene. 

(4)  IMPLICATION  FOR  FURTHER  RESEARCH 

The  tools  developed  facilitate  simulations  far  from  equilibrium.  The  result 
is  enhanced  correlative  and  predictive  capabilities  for  systems  undergoing  nano¬ 
second  to  picosecond  thermal  and  mechanical  processes.  This  is  a  rapidly-changing 
field  fundamental  to  the  understanding  of  complex  material  behavior. 
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I.  INTRODUCTION 


Systems  far  from  equilibrium  can  only  be  simulated  with  dynamic  equations  of 
motion  incorporati ng  nonequilibrium  temperature,  velocity,  and  stress  distribu¬ 
tions.  During  the  past  decade  a  variety  of  nonequilibrium  simulation  techniques 
have  been  developed.  During  the  past  three  years  a  consolidation  and  simplifica¬ 
tion  has  taken  place:  it  has  been  realized  that  simple  feedback  analysis  can  be 
applied  to  hydrodynamic  and  solid  mechanics  problems  in  a  way  consistent  with 
macroscopic  constitutive  analyses.  This  field  has  been  reviewed  recently^*2.3< 
The  highlights  are  as  follows: 

(1)  There  is  an  old  mechanical  principle,  due  to  Gauss,  "Gauss'  Principle  of 
Least  Constraint",  which  is  more  general  than  Newtonian  or  Lagrangian 
mechanics.  Gauss'  Principle  makes  it  possible  to  develop  equations  of 
motion  with  nonholonomic  constraints.  Such  constraints  include,  for 
instance,  the  restrictions  that  a  particular  process  be  carried  out  with 
a  given  time  and  space  variation  of  temperature  or  energy  or  stress. 
Gauss'  Principle  has  been  widely  tested,  during  the  past  three  years, 
and  found  to  be  useful  in  studying  steady  flow  far  from  equilibrium. 

(2)  There  is  a  new  idea,  due  to  Shuichi  Nosfe( Yokahama) ,  which  makes  it  possi 
ble  to  simulate,  with  molecular  dynamics,  systems  described  with  indepen 
dent  variables  such  as  temperature  and  stress,  as  opposed  to  energy  and 
volume. Nose's  ideas  are  extremely  interesting  for  theoretical 
mechanics  and  suggest  new  practical  techniques,  different  to  Gauss',  for 
treating  systems  far  from  equilibrium. 

(3)  Both  Gauss'  and  Nosfe's  ideas  can  be  viewed  as  special  cases  of  Control 
Theory^.  This  realization  has  increased  the  multiplicity  of  algorithms 
available,  but  with  the  advantage  of  a  systematic  description.  The 


consequences  of  this  development  are  now  only  beginning  to  be  explored. 
This  is  a  most  exciting  development  in  mechanics,  with  implications  for 
quantum  as  well  as  classical  systems. 

During  the  past  three  years  we  have  contributed  to  the  development  of  these 
ideas  while  carryi ng  out  an  investigation  of  solid  hexani trobenzene.  The  complex¬ 
ities  associated  with  modelling  this  material  precluded  a  complete  study  of  its 
response  to  shockwave  initiation  and  detonation,  but  the  groundwork  and  the  tools 
necessary  to  this  simulation  are  now  in  place. 

The  following  sections  describe  the  new  techniques  available  for  atomistic 
simulations,  a  simple  model  for  crystalline  hexani trobenzene,  and  recommendations 
for  further  work. 
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II.  NEW  TECHNIQUES  IN  SIMULATIONS  FAR  FROM  EQUILIBRIUM 


In  1972  Ashurst  simulated  the  flow  of  momentum  p  and  the  flow  of  heat  between 
atomistic  reservoirs.  He  performed  a  momentum  rescaling 

p’  =  (1  +  e2)  p  +  e1 

in  which  the  small  numbers  e-|  and  e2  were  chosen  to  restore  the  first  and 
second  moments  of  the  velocity  distribution  to  desired  values.  His  were  the 
first  systematic  studies  of  nonequilibrium  systems  using  nonequilibrium  molecular 

Q 

dynamics  . 

In  1982  Evans  and  Hoover  noted  that  the  rescaling  of  momenta  could  be  carried 
out  in  a  continuous  way  by  adding  a  constraint  force  to  the  equation  of  motion  for 
each  particle: 


Fc  =  -  2P  *  F0  • 

Hoover  pointed  out  that  this  same  constraint  force  can  be  derived  from  Gauss  Prin¬ 
ciple  of  Least  Constraint.  By  1985  it  became  clear  that  many  forms  of  nonequilib¬ 
rium  constraint  forces  are  possible.  These  were  applied  to  a  variety  of  systems 

3 

far  from  equi 1 i bri urn  . 

The  theoretical  bases  of  these  constraint  forces  advanced  too.  It  was  pos¬ 
sible  to  show  that  exact  linear  transport  coefficients  could  be  generated  in 
steady  nonequilibrium  simulations.  Even  nonlinear  effects,  such  as  the  decrease 
of  viscosity  with  strain  rate  and  the  increase  or  decrease  of  conductivity  with 
heat  flux,  depended  only  weakly  on  the  particular  constraint  force  used.  In  the 
present  section  we  summarize  the  two  successful  approaches  which  emerged  from  this 
work  and  point  out  their  usefulness  in  studying  systems  far  from  equilibrium. 


A.  Gauss'  Principle 


1 

Gauss  stated  that  any  conceivable  constraint  on  a  dynamical  system  should  be 
implemented  by  minimizing  the  mean  squared  value  of  the  constraint  force  divided 
by  the  corresponding  mass: 

o 

ZFc/m  minimized.  (Gauss'  Principle) 

This  principle  provides  a  unique  prescription  for  carrying  out  simulations  of 
nonequilibrium  flows.  Several  consequences  of  this  principle  have  been  explored. 
Constant-energy  and  constant-temperature  (with  temperature  T  -  p  /3mk,  where  k 
is  Boltzmann's  constant)  shear  flows  yield  nearly  identical  viscosity  coefficients 

Q 

even  far  from  equilibrium.  Self-diffusion,  shear  and  bulk  viscosity,  and  heat 
conduction  can  be  studied  in  very  small  systems  (two  particles  for  the  first  three 
transport  coefficients,  and  three  for  heat  flow).10-12  When  periodic  bound¬ 
aries  are  used,  to  eliminate  surface  effects,  the  resulting  two-particle  transport 
coefficients  lie  within  a  factor  of  two  of  the  many-body  thermodynamic  limit. 

Thus  it  appears  that  nonequilibrium  results  have  an  intrinsic  number-dependence  no 
greater  than  the  equilibrium  dependence,  1/N,  where  N  is  the  number  of  particles. * 

Gauss'  principle  can  lead  to  erroneous  results.  If  a  simulation  is  carried 
out  with  fixed  heat-flux  vector  and  fixed  temperature  the  corresponding  thermal 
conductivity  can  be  in  error,  in  a  dense  fluid,  by  several  percent. 

Nevertheless,  a  variety  of  simulations  led  to  accurate  values  for  the  trans¬ 
port  coefficients  of  simple  fluids  and  solids.  These  new  data  suggested  useful 
corresponding  states  relations  between  the  transport  coefficients  and  equilibrium 
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These  same  calculations  also  suggest  that  the  nonequilibrium  conditions  pre¬ 
vailing  in  a  strong  shockwave  be  simulated  by  dynamical  adiabatic  compressions, 
with  the  time  of  the  compression  chosen  to  correspond  to  the  picosecond  timescale 
of  rising  stress  in  a  real  shockwave.  In  such  a  rapid  compression  it  can  be  an¬ 
ticipated  that  coupling  of  the  vibrational  and  translational  modes  is  incomplete 
and  that  highly  nonequilibrium  energy  distributions  resulting  from  the  lack  of 
coupling  produce  chemical  reactions  at  far  different  rates  than  would  be  predicted 
from  Arhenius  kinetics.  The  relaxation  process  is  just  beginning  to  be  studied  in 
a  serious  way  by  using  nonequilibrium  molecular  dynamics. 16  These  results 
indicate  that  Nose's  Dynamics  is  best  suited  to  simulations  far  from  equilibrium. 


B.  Nose's  Dynamics 

Nose  pointed  out  that  if  the  friction  coefficient  z  in  the  generalized 
equation  of  motion 


mr  =•  F(r)  -  zp  , 


relaxes  according  to  a  first-order  differential  equation 


z  -  a(p2-3mkT)  , 


where  a  is  fixed,  the  equilibrium  distribution  is  identical  to  Gibbs'  canonical 
distribution  for  the  temperature  T. 

This  approach  has  been  generalized  to  nonequilibrium  systems  and  appears  to 
be  very  successful.  Just  as  is  the  case  with  Gauss'  principle  it  is  possible  to 
analyze  systems  close  to  equilibrium  and  to  show  that  Nose's  dynamics  provides 


-32- 


exact  linear  transport  coefficients.  By  making  temperature  depend  upon  space  and 
time,  and  by  introducing  also  stress  as  an  independent  (tensor)  variable  it  is 
straightforward  to  generalize  Nose''s  approach  to  the  treatment  of  nonequilibrium 
problems  such  as  the  adiabatic  compression  of  hexani trobenzene.  These  simulations 
have  not  yet  been  carried  out,  because  they  are  relatively  complex,  but  the  tools 
for  doing  this  are  now  available.  The  model  for  crystalline  hexani trobenzene  is 


described  in  the  next  section. 


III.  HEXANITROBENZENE 

Hexanitrobenzene,  (CNC^)^,  is  a  crystalline  high  explosive,  about  twice 
the  density  of  water,  first  synthesized  and  characterized  in  Russia^. 

Initial  attempts  to  make  this  material  in  the  United  States  failed.  Later,  when 
it  became  clear  that  hexanitrobenzene  is  sensitive  both  to  water  and  to  light, 
military  interest  in  it  subsided.  Nevertheless,  from  the  scientific  standpoint 
this  substance  is  extremely  interesting.  Because  it  contains  no  hydrogen  it  is 
reasonable  to  use  a  classical  description  for  its  motion. 

The  molecules  have  sixfold  rotational  symmetry,  but  with  the  nitro  groups 
aligned  like  blades  of  a  propellor  to  alleviate  steric  repulsion  of  neighboring 
oxygens.  The  crystal  structure  was  described  in  1966  by  Akopyan,  Struchkov,  and 
Dashevskii.  For  reasons  that  are  not  clear  the  sixfold  symmetry  is  broken  in  the 
crystalline  phase.  The  molecules  are  arranged  in  layers  with  nearly  hexagonal 
symmetry,  with  each  near ly-equi 1 ateral  triangle  of  molecules  having  a  "long"  side 
about  one  percent  longer  than  its  two  "short"  sides.  These  nearly  hexagonal 
sheets  are  stacked,  in  a  nearly-face-centered-cubic  arrangement:  ABCABCABC. . . , 
with  a  spacing  between  planes  4.08A,  considerably  smaller  than  the  lengths  of  the 
triangle  legs,  9.05  and  9.13A. 

Discussions  with  quantum  chemists  indicated  that  a  reliable  simulation  of  the 
breakup  of  a  hexanitrobenzene  molecule  is  beyond  the  state  of  the  art.  A  likely 
mechanism  is  the  breakup,  due  to  bending,  of  one  of  the  six  nitro  groups.  Faced 
with  a  lack  of  detailed  information  on  the  breakup  mechanism  we  modelled  the 
molecule  as  six  carbon  atoms  and  six  nitro  groups,  with  6x3+6x3=36  degrees  of 
freedom  per  molecule.  Interactions  within  the  molecules  can  be  treated  in  the 
quasiharmonic  approximation.  Interactions  between  molecules  were  based  on  a 
truncated  Leonard- Jones  potential  to  which  a  Hooke' s-1 aw  potential  was  added  in 
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such  a  way  as  to  bring  both  the  potential  and  its  first  derivative  (minus  the 
force)  to  zero  at  the  potential  cutoff. 

Exploratory  calculations,  in  two  dimensions,  established  that  the  energy  was 
lowered  by  about  7 %  if  the  molecules  were  allowed  to  rotate.  In  these  calcula¬ 
tions  the  Lennard-Jones  potential  was  cut  off  at  twice  the  distance  to  the  poten¬ 
tial  minimum  which  was  in  turn  three  times  the  C-C  or  C-NC^  separation  distance. 

In  the  three-dimensional  case  both  the  nearly-hexagonal  and  nearly-face-centered 
forms  of  the  structure  were  studied.  The  least-energy  rotation,  just  over  8 
degrees,  was  approximately  equal  to  that  found  in  two  dimensions.  The  intermolec- 
ular  and  interplanar  spacings  were  8.75$  and  2.67$,  respectively.  The  correspond¬ 
ing  experimental  values  are  (9.05,  9.05,  9.13)$  and  4.08$.  The  energies  of  the 
two  forms  studied  agreed  to  four-figure  accuracy.  Because  the  range  of  the 
potential  exceeds  twice  the  interplanar  spacing  the  numbers  are  in  fact  slightly 
di fferent. 

A  final  investigation  was  designed  to  match  the  experimental  spacings  of 
hexani trobenzene,  9.08$  and  4.08A.  A  stress-free  unit  cell  was  constructed  using 
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a  potential  minimum  at  4.56A  and  a  C-C  +  C-N0£  separation  of  2.79A.  Again  the 
two  forms  of  crystal  agreed  in  energy  to  four-figure  accuracy.  A  well  depth  of 
9.6xl0-14  ergs  was  chosen  approximately  to  reproduce  the  estimated  heat  of 
vaporization  of  the  crystal.  A  detailed  report  of  these  calculations  is  in 

preparation '8. 


IV.  RECOMMENDATIONS 

It  is  recommended  that  a  study  of  the  normal-mode  energies,  bond  displace¬ 
ments,  and  angular  distortions  in  the  model  of  hexani trobenzene  be  studied  as  a 
function  of  initial  temperature,  shockwidth  and  strength  of  shockwave.  These 
calculations  would  be  useful  in  assessing  the  validity  of  competing  theories  of 
vibrational  relaxation,  as  described  by  Holian^  and  in  further  developing 
techniques  for  predicting  the  sensitivity  of  new  explosive  chemicals  to  detonation 
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